La solución más conocida para este problema es $x= e^{-\lambda}x_0$. Cuya gráfica es
import matplotlib.pyplot as plt
import numpy as np
plt.style.use('https://github.com/dhaitz/matplotlib-stylesheets/raw/master/pitayasmoothie-dark.mplstyle')
lamda = 5
x0 = 10
t = np.arange(0,10,0.01)
solrealy = np.exp(-lamda*t)*x0
plt.plot(t,solrealy)
[<matplotlib.lines.Line2D at 0x1876e5a8d50>]
def EF_caso1(lamda,h,N,x0):
return ((1-lamda*h)**N)*x0
def EB_caso1(lamda,h,N,x0):
return (x0*(1/(1+lamda*h))**N)
def CN_caso1(lamda,h,N,x0):
return (x0*((2-lamda*h)/(2+lamda*h))**N)
hs = [0.1,0.2,0.03,0.4,0.5,0.75,1]
ts = [np.arange(0,10,i) for i in hs]
lamda = 5
x0 = 10
mat_parte1 = np.ones([3,len(hs)]) # 3 metodos* n hs * n tiempo
mat_parte1 = [list(mat_parte1[i]) for i in range(3)]
metodos_num_caso1 = [EF_caso1,EB_caso1,CN_caso1]
metodos = ['EF','EB','CN']
num = 0
plt.figure()
fig,axs = plt.subplots(1,3,layout = 'constrained')
for metodo,ax in zip(metodos_num_caso1,axs.flat):
for aux in range(len(hs)):
h = np.abs(ts[aux][0]-ts[aux][1])
N = np.arange(0,len(ts[aux]))
mat_parte1[num][aux] = metodo(lamda,h,N,x0)
ax.plot(ts[aux],metodo(lamda,h,N,x0))
#ax.set_ylim([-10,10])
#Faltan los labels y titulos pero ya
ax.set_title(f'{metodos[num]}')
num +=1
#Gráficacion
leg = [f'h = {pp}' for pp in hs]
axs[2].legend(leg,loc='center left', bbox_to_anchor=(1, 0.5))
<matplotlib.legend.Legend at 0x18771927b50>
<Figure size 640x480 with 0 Axes>
En este caso, el análisis de estabilidad nos indiga que h's van a converger a la función. Aquí se puede ver cómo principlamente en Euler forward la estabilidad dada por cada h influye mucho en que tan cercana se vuelve mi aproximación a la solución real
ks = [1,2,5,10]
hs = np.arange(0+0.1,1,0.1)
x0 = 10
def EF_caso2(h,k,xn):
return xn + h*(-k*xn**2)
def EB_caso2(h,k,xn):
a = h*k
b = 1
c= -xn
return (-b+np.sqrt(b**2-4*a*c))/ (2*a)
def CN_caso2(h,k,xn):
a = (h/2)*k
b = 1
c= -k*xn+k*(h/2)*xn
return (-b+np.sqrt(b**2-4*a*c)) / (2*a)
ts = [np.arange(0,10,i) for i in hs]
mat_parte1 = np.ones([3,len(hs)]) # 3 metodos* n hs * k * n tiempo
mat_parte1 = [list(mat_parte1[i]) for i in range(3)]
metodos_num_caso1 = [EF_caso2,EB_caso2,CN_caso2]
plt.figure()
fig,axs = plt.subplots(4,3,layout = 'constrained',figsize = (10,10))
nombre_metodo = ['EF','EB','CN']
num2 = 0
for k in ks:
num = 0
for metodo in metodos_num_caso1:
for aux in range(len(hs)):
h = np.abs(ts[aux][0]-ts[aux][1])
l_aux = [0]*len(ts[aux])
#N = np.arange(0,len(ts[aux]))
for i in range(len(ts[aux])):
if i == 0:
l_aux[0] = x0
else:
l_aux[i] = metodo(h,k,l_aux[i-1])
mat_parte1[num][aux] = l_aux
# mat_parte1[num][aux] = metodo(h,1,mat_parte1[num][aux])
axs[num2][num].plot(ts[aux],l_aux,label = f'h = {hs[aux]}')
#Faltan los labels y titulos pero ya
axs[num2][num].set_title(f'Metodo = {nombre_metodo[num]} , k = {k}')
axs[num2][num].set_ylim([-5,10])
if num2 == 0 and num == 2:
axs[num2][num].legend(loc='center left', bbox_to_anchor=(1, 0.5))
num +=1
num2 += 1
C:\Users\alanu\AppData\Local\Temp\ipykernel_27132\2543124031.py:5: RuntimeWarning: overflow encountered in double_scalars return xn + h*(-k*xn**2)
<Figure size 640x480 with 0 Axes>
F = [[0,1],[-1,0]]
x0 = [2,0]
hs = [0.1,0.5,1]
ts = [np.arange(0,10,i) for i in hs]
ns = [len(t) for t in ts]
l_x1 = [[[0]*n for n in ns] for _ in range(2)]
l_x2 = [[[0]*n for n in ns] for _ in range(2)]
def EF_caso3(F,x0,h,n):
eigenvalues, eigenvectors = np.linalg.eig(F)
E = eigenvectors
D = np.diag(eigenvalues)
E_inv = np.linalg.inv(eigenvectors)
#A_diag = eigenvectors @ D @ np.linalg.inv(eigenvectors)
I = np.identity(2)
Xn = np.matmul( np.matmul(np.matmul(E,( np.linalg.matrix_power((I+h*D),n))),E_inv),x0)
return Xn
def EB_caso3(F,x0,h,n):
eigenvalues, eigenvectors = np.linalg.eig(F)
E = eigenvectors
D = np.diag(eigenvalues)
E_inv = np.linalg.inv(eigenvectors)
#A_diag = eigenvectors @ D @ np.linalg.inv(eigenvectors)
I = np.identity(2)
Xn = np.matmul( np.matmul(np.matmul(E, np.linalg.matrix_power(np.linalg.inv(I+h*D), n)) ,E_inv) ,x0)
return Xn
#Solución análitica
def sol_analitica(t):
return np.array([[2 * np.cos(t)], [2 * np.sin(t)]])
t_solr = np.arange(0,10,0.1)
aux_mat = sol_analitica(t_solr)
x1_analitico = aux_mat[0][0]
x2_analitico = aux_mat[1][0]
# Comparadas con la sol análitica
metodos = ['EF','EB']
plt.figure()
fig,axs = plt.subplots(3,2,layout = 'constrained',sharex= True)
num = 0
for metodo in [EF_caso3,EB_caso3]:
num2 = 0
for aux in range(len(hs)):
l_aux = [0]*ns[aux]
#Paso
for n in range(ns[aux]):
l_aux[n] = metodo(F,x0,hs[aux],n)
x1_num = np.real(np.array(l_aux).T[0])
x2_num = np.real(np.array(l_aux).T[1])
l_x1[num][num2] = x1_num
l_x2[num][num2] = x2_num
axs[num2][0].plot(ts[aux],x1_num,label = f'h = {hs[aux]} ,m = {metodos[num]}')
axs[num2][1].plot(ts[aux],x2_num)
axs[num2][0].plot(t_solr,x1_analitico,linewidth = 1.5,linestyle = 'dotted')
axs[num2][1].plot(t_solr,x2_analitico,linewidth = 1.5,linestyle = 'dotted')
axs[num2][0].set_title(f'X1 h = {hs[aux]}')
axs[num2][1].set_title(f'X2 h = {hs[aux]}')
num2 +=1
num+=1
axs[0][1].legend(['Forward','Análitica','Backward'])
<matplotlib.legend.Legend at 0x1877244a890>
<Figure size 640x480 with 0 Axes>
plt.figure()
fig,axs = plt.subplots(1,2,layout ='constrained',figsize = (8,5))
num = 0
for ax in axs.flat:
for j in range(3):
ax.plot(l_x1[num][j],l_x2[num][j])
num +=1
axs[0].plot(x1_analitico,x2_analitico,linewidth = 1.5,linestyle = 'dotted')
axs[1].plot(x1_analitico,x2_analitico,linewidth = 1.5,linestyle = 'dotted')
axs[0].set_title('Forward')
axs[1].set_title('Backward')
axs[1].legend(['h = 0.1','h = 0.5','h = 1','Análitico'])
axs[0].set_xlim([-5,5])
axs[0].set_ylim([-5,5])
plt.text(-0.4, 1.1, 'Espacio de Fase', horizontalalignment='center',
verticalalignment='center', transform=ax.transAxes,fontsize = 13)
Text(-0.4, 1.1, 'Espacio de Fase')
<Figure size 640x480 with 0 Axes>
import numpy as np
import matplotlib.pyplot as plt
# Parámetros y configuración inicial
h = 0.1
r = 1
Ng = 15
tiempo = np.arange(0, 10, h)
espacio = np.arange(0, Ng, r)
mat_parte3 = np.ones([len(tiempo), Ng])
def fun_inicial(x):
s_asterisco = 5
return (1 / np.cosh((x - s_asterisco) / 2)) ** 2
mat_parte3[0] = fun_inicial(espacio)
# Cálculo de la matriz
for j in range(len(tiempo) - 1):
for i in range(Ng):
i_p3 = (i + 3) % Ng
i_p1 = (i + 1) % Ng
i_m1 = (i - 1) % Ng
i_m3 = (i - 3) % Ng
if j == 0 and i == 0:
print(i_p3)
print(i_p1)
print(i_m1)
print(i_m3)
mat_parte3[j + 1][i] = mat_parte3[j][i] + (h / (2 * r)) * (mat_parte3[j][i_p3] - mat_parte3[j][i_m1]) + (h / (8 * r**3)) * (mat_parte3[j][i_p3] - 3 * mat_parte3[j][i_p1] + 3 * mat_parte3[j][i_m1] - mat_parte3[j][i_m3])
# Graficar la matriz
fig, ax = plt.subplots()
cax = ax.imshow(mat_parte3, extent=[espacio[0], espacio[-1],tiempo[0], tiempo[-1]], cmap='gray',vmin = 0,vmax = 3) #, origin='lower'
cbar = fig.colorbar(cax, ax=ax)
ax.set_xlabel('Espacio')
ax.set_ylabel('Tiempo')
ax.set_title('Solitron')
plt.show()
3 1 14 12
from matplotlib.animation import FuncAnimation
x = range(Ng)
y = mat_parte3
# Inicializar la figura y los ejes
fig, ax = plt.subplots()
ax.set_xlim(0, 15)
ax.set_ylim(0, 1)
# Inicializar la línea que se actualizará en la animación
line, = ax.plot([], [], lw=2)
# Función de inicialización
def init():
line.set_data([], [])
return line,
# Función de actualización de la animación
def update(frame):
x_data = range(Ng)
y_data = mat_parte3[frame]
line.set_data(x_data, y_data)
return line,
# Crear la animación
ani = FuncAnimation(fig, update, frames=len(tiempo), init_func=init, blit=True)
ani.save( 's.gif', writer='imagemagick')
plt.show()
MovieWriter imagemagick unavailable; using Pillow instead.